Reconstruction of a surface topography

ABSTRACT

A system and method is provided for reconstructing a surface of an object. A system includes an input for receiving a 2-dimensional grid of measurements representing a surface of an object. Each grid point of the 2-dimensional grid of measurements includes information on the slope of the surface in the x-direction and y-direction. A processor selects a 2-dimensional part of the grid and fits a corresponding part of the surface to the measurements of all grid points in the selected part. For each grid point of the selected part, the fitting is based on both the corresponding first and second slope information. An output of the system provides a representation of at least the reconstructed surface part to a user.

The invention relates to a method of reconstructing a surface of anobject, where the object is represented by a 2-dimensional grid ofmeasurements. For each grid point the measurements include correspondinginformation on a first slope of the surface in a first direction and asecond slope of the surface in a different second direction. Theinvention further relates to software for executing the method and to asystem for reconstructing the surface of an object.

Surfaces of objects, such as silicon wafers, can be reconstructed frommeasured slopes. The slopes are typically measured by scanning thesurface along lines and measuring the slope. The measurements result ina 2-dimensional grid, where for each grid point a slope in a‘horizontal’ direction and a slope in a ‘vertical’ direction is given.The measurements may be obtained using deflectometry, where light (e.g.from a laser) is projected onto the surface and an angle of reflectionis measured, providing information on the slope. FIG. 1 shows a typicalequidistant measurement grid with measurements points along therectangular coordinate lines.

For many applications the actual surface (i.e. the topography of thesurface) needs to be determined based on the measured slopes. Currentreconstruction algorithms are based on carrying out a line integralalong a path through the measurement grid. For example, by performingsuch a line integral along each ‘horizontal’ path (ie. each pathparallel to one of the coordinate lines) the topography can bereconstructed. For each point along the path, the line integral uses theslope measured at the grid point in the direction of the path. Ingeneral the measured slopes contain measurement errors. Integrationalong the path accumulates all measurements errors in the slopes in thedirection of the path. In this way, the topography of the surfacedetermined at a grid point is therefore not determined uniquely butdepends on the path chosen though the grid. For example, a line integralalong a closed path like the path from (i,j) to (i+1, j), to (i+1, j+1)to (i, j+1) and back to (i,j) will in general not give the startingvalue for the surface topography.

In general, the measurements will not result in the regular equidistantmeasurement grid of FIG. 1 with measurements points along therectangular coordinate lines. In practice, measurement scan lines arenot always straight but may be curved. If a measurement system is usedthat obtains slope information only one direction during a scan,scanning needs to be performed in both directions. In such a case,measurement points in one direction may not be perfectly aligned withthe measurement points in the other direction. To compensate forirregularities in the measurements, usually a calibration procedure isperformed to correct the position of the measurement points to obtainthe perfect equidistant rectangular measurement grid. Such calibrationis usually based on interpolation that in itself introduces furthererrors in the measurements and is time-consuming.

It is an object of the invention to provide an improved method ofreconstruction of a topography of a surface. It is a further object toprovide software for performing such a method and a system for executingthe method.

To meet the object of the invention, the method of reconstructing asurface of an object, where the object is represented by a 2-dimensionalgrid of measurements, where for each grid point the measurements includecorresponding information on a first slope of the surface in a firstdirection and a second slope of the surface in a different seconddirection, includes selecting a 2-dimensional part of the grid andfitting a corresponding part of the surface to the measurements of allgrid points in the selected part, where the fitting for each grid pointof the selected part is based on both the corresponding first and secondslope information. By performing a fitting based on all measured slopesof the selected part, the effect of a localized measurement error issignificantly reduced to mainly the area of the error. According to themethod, much more measurements are fully used for the reconstruction.Instead of performing an integration operation along a one-dimensionalpath, now a 2-dimensional fitting operation is performed, involving muchmore measured data. Moreover, instead of using only one measured slopeof each involved grid point, now both slopes are used. This enablessignificant reduction in propagation of measurement errors. Preferably,a user can indicate the area in which accurate reconstruction accordingto the invention is desired. Outside the selected area, a conventionalreconstruction may be carried out In a preferred embodiment, asdescribed in the dependent claim 4, the fitting according to theinvention is performed over substantially the entire surface representedby the grid measurements.

According to dependent claim 2, the fitting is performed though aleast-square minimization operation. This is an effective way ofminimizing fitting errors.

According to dependent claim 3, the least square minimization operationis performed by solving an equation that describes a shape of a soapfilm loaded with a pressure field equal to a divergence of a slopevector including the first and second slope information. The inventorhad the insight that in this way the surface fitting according to theinvention can be expressed in a way similar to describing a shape of asoap film loaded with a pressure field. This enables using known methodsfor determining such a soap film shape to determine the topography ofthe surface.

According to dependent claim 5, for each point of the grid the first andsecond slope are measured using deflectometry.

To meet the object of the invention, a computer program productoperative to cause a processor to perform the steps of the method asclaimed in claim 1.

To meet the object of the invention, a system for reconstructing asurface of an object includes an input for receiving a 2-dimensionalgrid of measurements representing a surface of an object, where for eachgrid point the measurements include corresponding information on a firstslope of the surface in a first direction and a second slope of thesurface in a different second direction; a processor for, under controlof a program, selecting a 2-dimensional part of the grid and fitting acorresponding part of the surface to the measurements of all grid pointsin the selected part, where the fitting for each grid point of theselected part is based on both the corresponding first and second slopeinformation; and an output for providing a representation of at leastthe reconstructed surface part.

According to dependent claim 8, the system includes a measurement unitfor measuring for each measurement point of a measurement grid thecorresponding first and second slope information.

According to dependent claim 9, the measuring is performed alongnon-straight lines; the measurement grid being directly used for thereconstruction. The method and system according to the invention performa fitting that works as long as a connectivity between measurementpoints (analogous to the nodes belonging to a finite element in a FEMmesh) can be established. It is not required that the grid is neatlyarranged, e.g. that the measurement points in the two directions form anequidistant x, y or r, .phi. grid). No calibration of the measurementgrid to a grid used for the reconstruction is required, avoiding theintroduction of additional errors.

According to dependent claim 10, the measurement unit includes adeflectometry measurement unit.

These and other aspects of the invention are apparent from and will beelucidated with reference to the embodiments described hereinafter.

In the drawings:

FIG. 1 shows a quadrilateral finite element mesh corresponding to themeasurement grid;

FIG. 2 shows a block diagram of a system in which the invention can beemployed;

FIG. 3 shows a finite element mesh consisting of triangular elements;

FIG. 4 illustrates a piecewise bilinear base function on element meshshown on the 4 elements;

FIGS. 5A and B shown an exemplary initial surface and reconstructedsurface obtained from exact slopes with the method according to theinvention;

FIG. 6 shows the discretization error, i.e. difference betweenreconstructed and exact prescribed surface; and

FIGS. 7A-7F show the exact slopes, slopes containing noise used toreconstruct the surface, reconstructed surface and error inreconstructed surface.

FIG. 2 shows a block diagram of a system in which the invention can beemployed. In this exemplary system, the reconstruction of the surfacetopography is performed by a processor 100 that is loaded with asuitable program. The program may be stored in a permanent storage 150,such as a hard disk, and be executed from random access memory, such asDRAM (not shown). The processor operates on measurements receivedthrough an input 150. The input may be received from a measuring system130. The measuring system 130 may be external to the system in whichcase the measurements may be received through a suitable communicationsystem (e.g. wired or wireless LAN or wide area network, such asInternet). Preferably, the measuring system is part of the system forreconstruction of the surface. Advantageously, the measuring system is adeflectometry system, for example like the one described in U.S. Pat.5,831,738, hereby included by reference. For the description it isassumed that measurements are provided for a 2-dimensional grid ofmeasurements points, where for each point a slope in the x-direction andthe y-direction is given. The method works equally well for othersuitable coordinate systems, e.g. using an r, φ grid It is not requiredthat the grid is neatly arranged, e.g. that the measurement points inthe two directions form an equidistant x, y grid; it is sufficient thatconnectivity can be established between measurement points (analogous tothe nodes belonging to a finite element in a FEM mesh).

The input 120 is also able to receive input from a user. Shown are inputdevices, such as a mouse 140 and keyboard 150. Output to a user may besupplied via an output device 160, that may include a display. Inparticular the input may enable the user to indicate an area of thesurface (and/or grid) in which accurate reconstruction of the surface isrequired. For other areas, conventional reconstruction may take place.Such conventional reconstruction may be used to reconstruct the edge ofthe selected area. In the remainder, the reconstruction method accordingto the invention will be described.

Reconstruction method:

Consider a domain D⊂R² on which the topography of an unknown surface ismeasured with a deflectometry measurement system. Using this measurementsystem, the slopes {right arrow over (N)}({right arrow over (x)}) of theunknown surface z=z({right arrow over (x)}) are measured. If nomeasurement errors would occur, the unknown function that describes thesurface satisfies{right arrow over (∇)}z({right arrow over (x)})={right arrow over(N)}({right arrow over (x)}), {right arrow over (x)}∈D  (1)Hence, the surface z({right arrow over (x)}) is determined apart from aconstant function z=z₀. A unique surface is obtained by prescribing thesurface at some point {right arrow over (x)}₀∈D, hence z({right arrowover (x)}₀)=z₀.

In general the measured slopes contain measurement errors. Then, it ismost likely that no function z({right arrow over (x)}) exists thatsatisfies(1), because in all cases where rot({right arrow over(N)})={right arrow over (∇)}×{right arrow over (N)}≠0 no function existsthat satisfies (1). However, it is possible to integrate (1) along apath Γ({right arrow over (x)},{right arrow over (x)}₀) from {right arrowover (x)}₀ to {right arrow over (x)} in the domain D to obtain valuesfor the surface topography

$\begin{matrix}{{z( \overset{arrow}{x} )} = {z_{0} + {\int_{\Gamma{({\overset{arrow}{x},{\overset{arrow}{x}}_{0}})}}{\overset{arrow}{N} \cdot \ {\mathbb{d}\overset{arrow}{s}}}}}} & (2)\end{matrix}$Note, that if rot({right arrow over (N)})≠0 these values are dependenton the path Γ({right arrow over (x)},{right arrow over (x)}₀) and aretherefore not unique. This is seen by carrying out the line integration(2) along a closed path Γ({right arrow over (x)}₀,{right arrow over(x)}₀) enclosing a surface Ω_(Γ) and using Stokes theorem hence

$\begin{matrix}{{\int_{\Gamma{({{\overset{arrow}{x}}_{0},\overset{arrow}{x}})}}{\overset{arrow}{N}\ {\mathbb{d}\overset{arrow}{s}}}} = {\int_{\Omega_{\Gamma}}{\int{{{{rot}( \overset{arrow}{N} )}\  \cdot {\overset{arrow}{n}}_{\Omega}}{\mathbb{d}a}}}}} & (3)\end{matrix}$where {right arrow over (n)}_(Ω) is the unit normal on Ω_(Γ) in R³.Next, a method is proposed to find a unique function that approximatelysatisfies (1). Hence, the reconstructed surface obtained by the proposedmethod is independent of the fact that rot({right arrow over (N)}) iszero or not.

The unknown surface topography is approximated by a function z({rightarrow over (x)}) that satisfies (1) in a least squares sense, hence itis the function with z({right arrow over (x)}₀)=z₀ that minimizes thefunctional

$\begin{matrix}{{J(z)} = {\int_{D}{{( {{\overset{arrow}{\nabla}z} - \overset{arrow}{N}} ) \cdot ( {{\overset{arrow}{\nabla}z} - \overset{arrow}{N}} )}{\mathbb{d}a}}}} & (4)\end{matrix}$where the integration is carried out over the measurement domain D. Inorder to be a minimum for every disturbance δz({right arrow over (x)})of z({right arrow over (x)}) the minimizing function z({right arrow over(x)}) has to satisfy

$\begin{matrix}{{{\delta\;{J(z)}} = {{\int_{D}{\{ {{( {{\overset{arrow}{\nabla}z} - \overset{arrow}{N}} ) \cdot {\overset{arrow}{\nabla}\delta}}\; z} \}\ {\mathbb{d}a}}} = 0}},{{z( {\overset{arrow}{x}}_{0} )} = z_{0}}} & (5)\end{matrix}$This last formulation of the solution of the minimization problem (4) isthe weak formulation of{right arrow over (∇)}·{right arrow over (∇)}z={right arrow over(∇)}·{right arrow over (N)}, {right arrow over (x)}∈D  (6)with boundary conditionz({right arrow over (x)} ₀)=z ₀ ; {right arrow over (n)}·{right arrowover (∇)}z={right arrow over (n)}·{right arrow over (N)}, {right arrowover (x)}∈∂D  (7)where {right arrow over (n)} is the unit normal at the boundary ∂D.

The formulation (6) is also directly found by taking the divergence ofthe left and right hand side of (1). The equation (6) also describes thesolution of the mechanical problem of determining the shape of a soapfilm loaded with a pressure distribution equal to {right arrow over(∇)}·{right arrow over (N)}. Note, that in the weak formulation (5) theorder of differentiation is one, while it is two in the formulation (6).Hence, solving (5) instead of (6) means that it is sufficient that theunknown function z({right arrow over (x)}) lies in the Sobolev space

$\begin{matrix}{{H^{1}(D)} = \{ {{z: Darrow\Re }❘{{\int_{D}{{{\overset{arrow}{\nabla}z} \cdot {\overset{arrow}{\nabla}\ z}}{\mathbb{d}a}}} < \infty}} \}} & \;\end{matrix}$and is not necessarily in C¹(D) (first order derivatives exist and arecontinuous). Furthermore, in the formulation (6) the first orderderivatives of the measured slope {right arrow over (N)} must exist.

The advantage of solving (5) instead of carrying out the pathintegration (2) is that:

-   -   Both measured slope signals are used simultaneously to determine        the topography of the measured surface.    -   All measured slope data is used to reconstruct the topography of        the measured surface.    -   The determined topography is independent of an integration path        which has to be chosen using the simple line integration        methods.    -   Analogous to the loaded soap film problem, the effect of errors        (inaccuracies) in the measured slope data on the determined        topography is mainly limited to the location of the inaccuracy.        With the line integration methods a local error propagates along        the remaining integration path and hence has a larger effect        (larger error) on the reconstructed topography.

Preferably, equation (5) is solved using a Finite Element Method (FEM).Of course other methods may be used to solve that equation, e.g. using aspectral method. Alternatively, equation (6) can be solved with e.g. afinite difference method. The formulation of the Finite Element Methodto solve (5) will be carried out in Cartesian coordinates. For othercoordinate systems the formulation is analogous to the presentlydiscussed one.

Consider a Cartesian coordinate system (x, y, z). Suppose that themeasurement grid D_(M)⊂D⊂R² consists of n·m discrete points (x_(ij),y_(ji)), i=1 . . . n, j=1 . . . m, with x_(i,j)<x_(i+1,j) andy_(j,i)<y_(j+1,i). With these discrete measurement points aquadrilateral finite element mesh is established by connecting thepoints (x_(i,j), y_(j,i)), (x_(i+1,j),y_(j,i+1)),(x_(i+1,j+1),y_(j+1,j+1)) and (x_(i,j+1), y_(i+1,j)) for i=1 . . . n-1and j=1 . . . m-1, as shown in FIG. 1. Hence, the measurement domain isdivided in (n-1)·(m-1) quadrilateral elements D_(e) using the nodes(x_(ij), y_(ji)) of the measurement grid. The choice of dividing themeasurement domain D into quadrilateral elements is not a restriction.For instance, if each quadrilateral element is divided in two triangles,as shown in FIG. 3, a finite element mesh consisting of triangularelements can be created. In that case the number of elements increaseswith a factor 2.

The unknown surface topography z=z(x, y) is written as a linearcombination of a finite set of approximating functions φ_(k)(x, y) k=1 .. . n·m, hence

$\begin{matrix}{{z( {x,y} )} = {\sum\limits_{k = 1}^{m \cdot n}{a_{k}{\varphi_{k}( {x,y} )}}}} & (8)\end{matrix}$Substitution of (8) in (3) and minimizing the functional J(z) withrespect to the coefficients a_(k), or equivalently substitution of (8)in (4) with δz=φ₁(x, y), yields

$\begin{matrix}{{{\sum\limits_{k = 1}^{m \cdot n}{\int_{D}{\int{{{\overset{arrow}{\nabla}\varphi_{k}} \cdot {\overset{arrow}{\nabla}\varphi_{l}}}{\mathbb{d}a}*a_{k}}}}} = {\int_{D}{\int{{\overset{arrow}{N} \cdot {\overset{arrow}{\nabla}\varphi_{i}}}\ {\mathbb{d}a}}}}},{l = {1\mspace{14mu}\ldots\mspace{14mu}{m \cdot n}}}} & (9)\end{matrix}$This singular set of linear equations is made regular by replacing thel^(th) equation, where l is the index at which |φ_(l)({right arrow over(x)}₀)| has a maximum, by the condition z({right arrow over (x)}₀)=z₀.The solution of this set of linear equations determines theapproximation (8) of the measured surface topography.

In the finite element approximation the subdivision of the measurementdomain D in the finite elements D_(e) is used to write the integrals in(9) as

$\begin{matrix}{{\int_{D}{\int{{{\overset{arrow}{\nabla}\varphi_{k}} \cdot {\overset{arrow}{\nabla}\varphi_{l}}}\ {\mathbb{d}a}}}} = {\sum\limits_{e}{\int_{D_{e}}{\int{{{\overset{arrow}{\nabla}\varphi_{k}} \cdot {\overset{arrow}{\nabla}\varphi_{l}}}{\mathbb{d}a}}}}}} & (10) \\{{\int_{D}{\int{{\overset{arrow}{N} \cdot {\overset{arrow}{\nabla}\varphi_{l}}}\ {\mathbb{d}a}}}} = {\sum\limits_{e}{\int_{D_{e}}{\int{{\overset{arrow}{N} \cdot {\overset{arrow}{\nabla}\varphi_{l}}}{\mathbb{d}a}}}}}} & (11)\end{matrix}$The functions φ_(k) are chosen to be bi-linear per element and have thevalue 1 in one node and 0 in all other nodes, as shown in FIG. 4. FIG. 4illustrates a piecewise bilinear base function on element mesh shown onthe 4 elements: [−1,0]×[−1,0], [0,1]×[−1,0], [0,1]×[0,1] and[−1,0]×[0,1].

For the evaluation of the integrals over the elements an iso-parametrictransformation of an actual quadrilateral element to the standardquadrilateral element is carried out

$\begin{matrix}{{\overset{arrow}{x}( \overset{arrow}{\xi} )} = {\sum\limits_{i = 1}^{4}{{\overset{arrow}{x}}_{k{(i)}}{\varphi_{k{(i)}}( \overset{arrow}{\xi} )}}}} & (11)\end{matrix}$where {right arrow over (x)}=(x,y)^(T), {right arrow over (ξ)}=(ξ,η)^(T)and the node number function k(i) is such that the actual 4 nodes of anelement are used. Then an integral of a function ƒ({right arrow over(x)}) over an element D_(e) is seen to be

$\begin{matrix}{{\int_{D_{e}}{\int{{f( \overset{arrow}{x} )}{\mathbb{d}a}}}} = {\int_{- 1}^{1}{\int_{- 1}^{1}{{f( {\overset{arrow}{x}( \overset{arrow}{\xi} )} )}{J( \overset{arrow}{\xi} )}{\mathbb{d}\xi}\;{\mathbb{d}\eta}}}}} & (12)\end{matrix}$where J({right arrow over (ξ)}) is the determinant of the Jacobianmatrix

$F = {\frac{\partial\overset{arrow}{x}}{\partial\overset{arrow}{\xi}}.}$The bilinear approximating or basis functions φ_(k)(ξ,η) are seen to be

$\begin{matrix}{{\varphi_{1}( {\xi,\eta} )} = {\frac{1}{4}( {1 - \xi} )( {1 - \eta} )}} & (13) \\{{\varphi_{2}( {\xi,\eta} )} = {\frac{1}{4}( {1 + \xi} )( {1 - \eta} )}} & (14) \\{{\varphi_{3}( {\xi,\eta} )} = {\frac{1}{4}( {1 + \xi} )( {1 + \eta} )}} & (15) \\{{\varphi_{4}( {\xi,\eta} )} = {\frac{1}{4}( {1 - \xi} )( {1 + \eta} )}} & (16)\end{matrix}$and their gradients are

$\begin{matrix}{{\frac{\partial{\varphi_{1}( {\xi,\eta} )}}{\partial\xi} = {\frac{- 1}{4}( {1 - \eta} )}};{\frac{\partial{\varphi_{1}( {\xi,\eta} )}}{\partial\eta} = {\frac{- 1}{4}( {1 - \xi} )}}} & (17) \\{{\frac{\partial{\varphi_{2}( {\xi,\eta} )}}{\partial\xi} = {\frac{1}{4}( {1 - \eta} )}};{\frac{\partial{\varphi_{2}( {\xi,\eta} )}}{\partial\eta} = {\frac{- 1}{4}( {1 + \xi} )}}} & (18) \\{{\frac{\partial{\varphi_{3}( {\xi,\eta} )}}{\partial\xi} = {\frac{1}{4}( {1 + \eta} )}};{\frac{\partial{\varphi_{3}( {\xi,\eta} )}}{\partial\eta} = {\frac{1}{4}( {1 + \xi} )}}} & (19) \\{{\frac{\partial{\varphi_{4}( {\xi,\eta} )}}{\partial\xi} = {\frac{- 1}{4}( {1 + \eta} )}};{\frac{\partial{\varphi_{4}( {\xi,\eta} )}}{\partial\eta} = {\frac{1}{4}( {1 - \xi} )}}} & (20)\end{matrix}$In the evaluation of the integrals in the right hand side of (9),integrals involving measured slope values occur. Using (10), (11) and(12) these integrals are transformed to integrals over the standardquadrilateral element. On the standard element the variation of theslope vector is approximated with a bilinear function

$\begin{matrix}{{\overset{arrow}{N}( {\overset{arrow}{x}( \overset{arrow}{\xi} )} )} = {\sum\;{{\overset{arrow}{N}( {\overset{arrow}{x}}_{k{(i)}} )}{\varphi_{k{(i)}}(\xi)}}}} & (21)\end{matrix}$Substitution of (21) in the expression obtained by (10), (11) and (12)yields integrals of the form

$\begin{matrix}{\int_{- 1}^{1}{\int_{- 1}^{1}{{\varphi_{k}( {\xi,\eta} )}\frac{\partial{\varphi_{l}( {\xi,\eta} )}}{\partial\xi}{J( {\xi,\eta} )}{\mathbb{d}\xi}\;{\mathbb{d}\eta}}}} & (22) \\{\int_{- 1}^{1}{\int_{- 1}^{1}{{\varphi_{k}( {\xi,\eta} )}\frac{\partial{\varphi_{l}( {\xi,\eta} )}}{\partial\eta}{J( {\xi,\eta} )}{\mathbb{d}\xi}\;{\mathbb{d}\eta}}}} & (23)\end{matrix}$These integrals can be evaluated analytically. However, in the presentwork a 4 point Gaussian quadrature rule is applied. The integrationpoints and weights of this method are

$\begin{matrix}{{( {\xi,\eta} )_{1} = ( {{\frac{- 1}{3}\sqrt{3}},{\frac{- 1}{3}\sqrt{3}}} )};{w_{1} = 1}} & (24) \\{{( {\xi,\eta} )_{2} = ( {{\frac{1}{3}\sqrt{3}},{\frac{- 1}{3}\sqrt{3}}} )};{w_{2} = 1}} & (25) \\{{( {\xi,\eta} )_{3} = ( {{\frac{1}{3}\sqrt{3}},{\frac{1}{3}\sqrt{3}}} )};{w_{3} = 1}} & (26) \\{{( {\xi,\eta} )_{4} = ( {{\frac{- 1}{3}\sqrt{3}},{\frac{1}{3}\sqrt{3}}} )};{w_{4} = 1}} & (27)\end{matrix}$

Carrying out all integrations, the element matrices and vectors areassembled in the large matrix and vector. Now the condition z({rightarrow over (x)}₀)=z₀ is incorporated into the equations, as discussedbefore and the resulting set of linear discrete equationsA_(lk)a_(k)=ƒ_(l)  (28)has to be solved. Several techniques exist to solve (28). It ispreferred to use a Gaussian elimination procedure (LU decomposition),which is a direct method. If the system of equations is very large,iterative procedures can be used to solve (28), e.g. Gauss-Seidel,Multigrid, etc.

EXAMPLES

In this section the accuracy and the efficiency of the proposed methodis investigated using virtual measurement data.

Suppose the measurement domain is D=[0,1]×[0,1]⊂R². The test surface isz(x, y)=sin²(2πx)sin²(2πy)  (29)with exact slope vector

$\begin{matrix}{{\overset{arrow}{N}\begin{pmatrix}{{\partial z}\text{/}{\partial x}} \\{{\partial z}\text{/}{\partial y}}\end{pmatrix}} = \begin{pmatrix}{4{{\pi sin}( {2\pi\; x} )}{\cos( {2\pi\; x} )}{\sin^{2}( {2\pi\; y} )}} \\{4{{\pi sin}^{2}( {2\pi\; x} )}{\sin( {2\pi\; y} )}{\cos( {2\pi\; y} )}}\end{pmatrix}} & (30)\end{matrix}$The initial surface and reconstructed surface obtained from the exactslopes (30) with the presently proposed method using an equidistantmeasurement grid of 21×21 points are shown in FIGS. 5A and B,respectively. The discretization error, i.e. difference betweenreconstructed and exact prescribed surface, is shown in FIG. 6. Thereconstructed surface approximates the initial surface very good. Next,the accuracy of the obtained reconstructed surface is investigated whenthe measurement grid is varied, and/or noise is added to the slopevector.

First, the discretization error to reconstruct the prescribed testsurface with a finite discrete mesh is investigated. The discretizationerror to reconstruct the prescribed test surface with a finite discretemesh are shown in Table 1, where ∥ƒ∥_(L) ₂ is the L₂ norm

$\begin{matrix}{{f}_{L_{2}} = \sqrt{\int{\int_{D}{f^{2}{\mathbb{d}a}}}}} & (31)\end{matrix}$and ∥f∥_(∞) is the infinity (or max) norm

$\begin{matrix}{{f}_{\infty} = {\max\limits_{{\overset{arrow}{x}\;}_{l,j}}{{f( {\overset{arrow}{x}}_{l,j} )}}}} & (32)\end{matrix}$

The discretization error reduces quadratically with the mesh size (i.e.a typical dimension of and element, e.g. the radius of the largestcircle that fits in an element). This behavior is expected using thebi-linear elements.

TABLE 1 Discretization error to reconstruct the prescribed test surface.n m ∥z_(h) − z_(exact)∥_(∞) ∥z_(h) − z_(exact)∥_(L) ₂ 10 10 76.808e-337.788e-3 20 20 24.838e-3  9.2565e-3 40 40  6.1787e-3  2.3027e-3 80 80 1.5428e-3  0.57495e-3

Secondly, using the same mesh (measurement grid) white noise (randomvalues) with an amplitude equal to the amplitude of the slope values isadded to the slope vector

$\begin{matrix}{{ {\overset{arrow}{N}( \overset{arrow}{x} )}\longrightarrow{\overset{arrow}{N}( \overset{arrow}{x} )}  + {{a\begin{pmatrix}{r_{x}( \overset{arrow}{x} )} \\{r_{y}( \overset{arrow}{x} )}\end{pmatrix}}\mspace{14mu}{where}\mspace{14mu} a}} = {{\max\limits_{\overset{arrow}{x},i}( {N_{i}( \overset{arrow}{x} )} )} - {\min\limits_{\overset{arrow}{x},i}( {N_{i}( \overset{arrow}{x} )} )}}} & (33)\end{matrix}$and r_(x)({right arrow over (x)}) and r_(y)({right arrow over (x)})generate random values in the range

$\lfloor {{- \frac{1}{2}},\frac{1}{2}} \rfloor.$In this case the difference between the reconstructed and the initialtest surface is caused by a discretization error and by the noise in theslope data. FIG. 7 shows: Original slopes: (A) in x-direction, (B) iny-direction; used slopes to reconstruct surface: (C) slope with noise inx-direction, (D) slope with noise in y-direction; and reconstructedsurface (E) and error in reconstructed surface (F), for-a gridconsisting of 21×21 points. Even using the slopes with the added whitenoise having the same amplitude as that of the original slopes, thereconstructed surface approximates the original surface fairly good. Asthe problem is linear in the known slopes, the difference caused by thenoise in the slope vector is linear with the amplitude of the noise. Thedifference between the reconstructed surface and the initial testsurface is shown in Table 2. As follows from Table 1 and Table 2, thedifference between the reconstructed and original surfaces is mainlycaused by the noise in the slopes. If the amplitude of the noise ishalved, also the difference between the reconstructed and originalsurface halves. Note, that the behavior of the noise in the twoconsidered cases of different amplitudes is different, because of thenot repeatable random generated values.

TABLE 2 Difference between reconstructed and initial test surface, whennoise is added to slope vectors. noise amplitude = a noise amplitude =0.5a n m ∥z_(h) − z_(exact)∥_(∞) ∥z_(h) − z_(exact) ∥_(L) ₂ ∥z_(h) −z_(exact)∥_(∞) ∥z_(h) − z_(exact)∥_(L) ₂ 10 10 0.69997 0.24197 0.284330.08341 20 20 0.41839 0.12651 0.216 0.061669 40 40 0.31338 0.0755380.13551 0.042426 80 80 0.18825 0.050282 0.084235 0.018521

It should be noted that the above-mentioned embodiments illustraterather than limit the invention, and that those skilled in the art willbe able to design many alternative embodiments without departing fromthe scope of the appended claims. In the claims, any reference signsplaced between parentheses shall not be construed as limiting the claim.The words “comprising” and “including” do not exclude the presence ofother elements or steps than those listed in a claim. The invention canbe implemented by means of hardware comprising several distinctelements, and by means of a suitably programmed computer. Where thesystem/device/apparatus claims enumerate several means, several of thesemeans can be embodied by one and the same item of hardware. The computerprogram product may be stored/distributed on a suitable medium, such asoptical storage, but may also be distributed in other forms, such asbeing distributed via the Internet or wireless telecommunicationsystems.

1. A method of reconstructing a surface of an object; the object beingrepresented by a 2-dimensional grid of measurements, where for each gridpoint the measurements include corresponding information on a firstslope of the surface in a first direction and a second slope of thesurface in a different second direction; the method including: in aprocessor, selecting a 2-dimensional part of the grid over which anaccurate reconstruction may be carried out, fitting a corresponding partof the surface to the measurements of all grid points in the selectedpart, thereby significantly reducing the effect of a localizedmeasurement error to the area of the selected 2-dimensional part, andproviding a representation of at least the reconstructed surface part,where the fitting for each grid point of the selected part is based onboth the corresponding first and second slope information.
 2. A methodas claimed in claim 1, including performing the fitting through aleast-square minimization operation.
 3. A method as claimed in claim 2,including performing the least square minimization operation by solvingan equation that describes a shape of a soap film loaded with a pressurefield equal to a divergence of a slope vector including the first andsecond slope information.
 4. A method as claimed in claim 1, wherein theselected part of the grid is substantially the entire grid.
 5. A methodas claimed in claim 1, including measuring for each point of the gridthe first and second slope using deflectometry.
 6. A computer-readablestorage medium comprising computer-readable instructions operative tocause a processor to perform the steps of the method as claimed inclaim
 1. 7. A system for reconstructing a surface of an objectincluding: an input for receiving a 2-dimensional grid of measurementsrepresenting a surface of an object, where for each grid point themeasurements include corresponding information on a first slope of thesurface in a first direction and a second slope of the surface in adifferent second direction; a processor, under control of a program, for(a) selecting a 2-dimensional part of the grid over which an accuratereconstruction may be carried out, and (b) fitting a corresponding partof the surface to the measurements of all grid points in the selectedpart, thereby significantly reducing the effect of a localizedmeasurement error to the area of the selected 2-dimensional part, wherethe fitting for each grid point of the selected part is based on boththe corresponding first and second slope information; and an output forproviding a representation of at least the reconstructed surface part.8. A system as claimed in claim 7, wherein the system includes ameasurement unit for measuring for each measurement point of ameasurement grid the corresponding first and second slope information.9. A system as claimed in claim 8, wherein the measuring is performedalong non-straight lines; the measurement grid being directly used forthe reconstruction.
 10. A system as claimed in claim 8, wherein thesystem the measurement unit includes a deflectometry measurement unit.11. A method of reconstructing a surface of an object; the object beingrepresented by a 2-dimensional grid of measurements, where for each gridpoint the measurements include corresponding information on a firstslope of the surface in a first direction and a second slope of thesurface in a different second direction; the method including: in aprocessor selecting a 2-dimensional part of the grid over which anaccurate reconstruction may be carried out, and fitting a correspondingpart of the surface to the measurements of all grid points in theselected part, thereby significantly reducing the effect of a localizedmeasurement error to the area of the selected 2-dimensional part, wherethe fitting for each grid point of the selected part is based on boththe corresponding first and second slope information, whereby saidfitting is performed through a least-square minimization operation bysolving an equation that describes a shape of a soap film loaded with apressure field equal to a divergence of a slope vector including thefirst and second slope information, and providing a representation of atleast the reconstructed surface part.
 12. A system for reconstructing asurface of an object including: an input for receiving a 2-dimensionalgrid of measurements representing a surface of an object, where for eachgrid point the measurements include corresponding information on a firstslope of the surface in a first direction and a second slope of thesurface in a different second direction; a processor, under control of aprogram, for selecting a 2-dimensional part of the grid over which anaccurate reconstruction may be carried out, and fitting a correspondingpart of the surface to the measurements of all grid points in theselected part, thereby significantly reducing the effect of a localizedmeasurement error to the area of the selected 2-dimensional part, wherethe fitting for each grid point of the selected part is based on boththe corresponding first and second slope information; and an output forproviding a representation of at least the reconstructed surface part,wherein the system includes a measurement unit for measuring for eachmeasurement point of a measurement grid the corresponding first andsecond slope information and wherein the measurement unit includes adeflectometry measurement unit.
 13. A computer program productcomprising one or more computer-readable storage media having thereoncomputer executable instructions that, when executed by one or moreprocessors of a system for reconstructing the surface of an object,cause the system to perform the following: select a 2-dimensional partof the grid over which an accurate reconstruction may be carried out;fit a corresponding part of the surface to the measurements of all gridpoints in the selected part, thereby significantly reducing the effectof a localized measurement error to the area of the selected2-dimensional part; and provide a representation of at least thereconstructed surface part, where the fitting for each grid point of theselected part is based on both the corresponding first and second slopeinformation.